{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Survival Analysis\n",
    "\n",
    "- The first function draws the Survival Curve (Kaplan-Meier curve).\n",
    "- The second function implements the logrank test, comparing two survival curves.\n",
    "\n",
    "The formulas and the example are taken from Altman, Chapter 13\n",
    "\n",
    "Author : Thomas Haslwanter, Date : Feb-2017"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "from urllib.request import urlopen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def kaplanmeier(data):\n",
    "    '''Determine and the Kaplan-Meier curve for the given data.\n",
    "    Censored times are indicated with \"1\" in the second column, uncensored with \"0\"'''\n",
    "    times = data[:,0]\n",
    "    censored = data[:,1]\n",
    "    atRisk = np.arange(len(times),0,-1)\n",
    "    \n",
    "    failures = times[censored==0]\n",
    "    num_failures = len(failures)\n",
    "    p = np.ones(num_failures+1)\n",
    "    r = np.zeros(num_failures+1)\n",
    "    se = np.zeros(num_failures+1)\n",
    "    \n",
    "    # Calculate the numbers-at-risk, the survival probability, and the standard error\n",
    "    for ii in range(num_failures):\n",
    "        if failures[ii] == failures[ii-1]:\n",
    "            r[ii+1] = r[ii]\n",
    "            p[ii+1] = p[ii]\n",
    "            se[ii+1] = se[ii]\n",
    "            \n",
    "        else:\n",
    "            r[ii+1] = np.max(atRisk[times==failures[ii]])\n",
    "            p[ii+1] = p[ii] * (r[ii+1] - sum(failures==failures[ii]))/r[ii+1]\n",
    "            se[ii+1] = p[ii+1]*np.sqrt((1-p[ii+1])/r[ii+1])\n",
    "            # confidence intervals could be calculated as ci = p +/- 1.96 se\n",
    "    \n",
    "    # Plot survival curve (Kaplan-Meier curve)\n",
    "    # Always start at t=0 and p=1, and make a line until the last measurement\n",
    "    t = np.hstack((0, failures, np.max(times)))\n",
    "    sp = np.hstack((p, p[-1]))\n",
    "    \n",
    "    return(p,atRisk,t,sp,se)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Survival Probability')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADt0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjByYzIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy/EUOrgAAAaYUlEQVR4nO3de5gcdZ3v8fcnk+AABhEkXDIggQ2XRJIRhuuKC7JnE0TNQWSF7Aa87JnlSFZZdtH4eFbJqo+ssoAaEHMgJqDcFM4SOQjmoFyezSJM3BAI4RJJgCHxJMRjwj237/mjakJn0t1TM+nqnp76vJ5nnu6q+lXN9weT/nb96ndRRGBmZsU1rNEBmJlZYzkRmJkVnBOBmVnBORGYmRWcE4GZWcE5EZiZFVxuiUDSHElrJD1R4bgkfU/ScklLJB2dVyxmZlZZnncEc4HJVY6fDoxNfzqBH+QYi5mZVZBbIoiIB4E/VCkyBbghEg8De0raP694zMysvOEN/N2jgRdLtrvTfat7F5TUSXLXwO67737MEUcc0e9f9tryFQx7841+nfOOjW/x1i7vYPW+7+337xsK3ti0hV1HtHDIPrs3OhQz20mLFi16OSL2KXeskYlAZfaVne8iImYDswE6Ojqiq6srz7jedsopyev999fn9w0yn/zhfwBw69+e2OBIzGxnSXq+0rFG9hrqBg4s2W4DVjUoFjOzwmpkIpgPnJf2HjoBWB8ROzQLmZlZvnJrGpJ0M3AK8B5J3cDXgBEAEXEtcDfwYWA58Drw6bxiMTOzynJLBBFxbh/HA7gwr99vZsWzadMmuru7efPNNxsdSsO0trbS1tbGiBEjMp/TyIfFZmY11d3dzciRIzn44IORyvVHGdoignXr1tHd3c2YMWMyn+dE0JfFi9/uPdRfU6dCZ2dNwzGzyt58883CJgEASey9996sXbu2X+c5EVQzderAz128OHl1IjCrq6ImgR4Dqb8TQTWdnQP/IB/oXYSZWZ05EVhVT67esG1gWTOb0j6aqccf1OgwrABaWlo46qij2LRpE8OHD+f888/noosuYtiwyr31V65cycKFC5maoRVi8uTJPPzww3zgAx/grrvuqknMnobaKprSPppx++/R6DB22pOrN3Dn4pcaHYYVxK677srixYtZunQpCxYs4O6772bmzJlVz1m5ciU33XRTputfcskl3HjjjbUIdRvfEVhFU48/aEh8ix4KdzTWnEaNGsXs2bM59thjufTSS3n++eeZNm0ar732GgCzZs3ipJNOYsaMGSxbtoz29nbOP/98zjzzzLLlAE477TTur/G0N04EZjYkzfz5Up5ctaGm1xx3wB587aPj+3XOIYccwtatW1mzZg2jRo1iwYIFtLa28uyzz3LuuefS1dXFZZddxuWXX76tqef1118vWy4vTgRmZjlLxs8mA96mT5/O4sWLaWlp4ZlnnilbPmu5WnEiMLMhqb/f3PPy3HPP0dLSwqhRo5g5cyb77rsvjz32GFu3bqW1tbXsOVdeeWWmcrXih8VmZjlZu3YtF1xwAdOnT0cS69evZ//992fYsGHceOONbNmyBYCRI0fyyiuvbDuvUrm8+I7AzKyG3njjDdrb27d1H502bRoXX3wxAJ/73Oc466yz+OlPf8qpp57K7rsniz5NmDCB4cOHM3HiRD71qU9VLAdw8skn89RTT/Hqq6/S1tbG9ddfz6RJk3YqZicCM7MaqvbtfezYsSxZsmTb9re+9S0ARowYwX333bdd2XLlAB566KFahbqNm4bMzArOdwRWCKUjpD3K2Gx7TgQ25E1pH73t/ZOrk37lTgRmb3MisCGvdIS0Rxmb7ciJIE/l1jLwGgVmNsg4EeSl3CyCXqPAzAYh9xrKS2cn3H//9j/t7Y2Nycxy19LSQnt7O+PHj2fixIlcccUVbN26teo5WWcfXbx4MSeeeCLjx49nwoQJ3HrrrTWJ2XcEZmY11DMNNcCaNWuYOnUq69evrzoVdU8i6Gs9gt12240bbriBsWPHsmrVKo455hgmTZrEnnvuuVMx+47AzCwnPdNQz5o1i4hg5cqVnHzyyRx99NEcffTRLFy4EIAZM2bw0EMP0d7ezpVXXlmx3GGHHcbYsWMBOOCAAxg1alS/1ycux3cEZjY0XXTR28/laqW9Ha66ql+n5DUN9SOPPMLGjRs59NBDd7paTgRmZjmr9TTUq1evZtq0acybN6/qEphZORGY2dDUz2/uean1NNQbNmzgjDPO4Bvf+AYnnHBCTWJ0Iqi3cmML8uRxCzsonW7CsvPUHP1XbhrqtrY2hg0bxrx586pOQ12u3MaNGznzzDM577zzOPvss2sWpxNBPfXRI6DmPG5hB6XTTVh2npojuzynob7tttt48MEHWbduHXPnzgVg7ty5tO9k13T1tF01i46Ojshz7c4hpefOo8YLXVvx9NxB3fq3JzY4kuqWLVvGkUce2egwGq7cfwdJiyKio1x5dx81Mys4JwIzs4JzIjCzIaXZmrtrbSD1dyIwsyGjtbWVdevWFTYZRATr1q2r2C21EvcaMrMho62tje7u7ppMu9CsWltbaWtr69c5TgRDXb3HLVhlHtORuxEjRjBmzJhGh9F0nAiGsnqPW7DKPKbDBrFcE4GkycB3gRbguoi4rNfxdwE/Bg5KY7k8In6UZ0yF0tnpD57BwndlNojl9rBYUgtwNXA6MA44V9K4XsUuBJ6MiInAKcC/Stolr5jMzGxHed4RHAcsj4jnACTdAkwBniwpE8BISQLeCfwB2JxjTGY2QPWao8lzGtVfnolgNPBiyXY3cHyvMrOA+cAqYCTwyYjYYU03SZ1AJ8BBB/kPxKze6jVHk+c0aow8E4HK7OvduXcSsBj4EHAosEDSQxGxYbuTImYDsyGZayiHWM2siqnHH1SXD2fPCtsYeQ4o6wYOLNluI/nmX+rTwB2RWA6sAI7IMSYzM+slzzuCR4GxksYALwHnAL37M74AnAY8JGlf4HDguRxjMmuccmM6PLbABoHcEkFEbJY0HbiXpPvonIhYKumC9Pi1wNeBuZIeJ2lK+lJEvJxXTGYNU25Mh8cW2CCR6ziCiLgbuLvXvmtL3q8C/iLPGMwGhXJjOjy2wAYJTzpnZlZwTgRmZgXnRGBmVnBOBGZmBefZR81sUKnXVBbNaNwBe/C1j46v+XWdCMxs0KjXVBa2PScCs0bKunBQQQae1WsqC9ueE4FZo2RdOMgDzyxnTgRmjZJ14SAPPLOcudeQmVnBORGYmRWcE4GZWcE5EZiZFVyfiUBSl6QLJb27HgGZmVl9ZbkjOAc4AHhU0i2SJqWLzZuZ2RDQZyKIiOUR8RXgMOAmYA7wgqSZkvbKO0AzM8tXpmcEkiYA/wp8B7gd+ASwAfhVfqGZmVk99DmgTNIi4I/A9cCMiHgrPfQbSX+aZ3BmZpa/LCOLz46I7RaUlzQmIlZExMdzisvMzOokS9PQzzLuMzOzJlTxjkDSEcB44F2SSr/57wG05h2YmZnVR7WmocOBjwB7Ah8t2f8K8N/yDMrMzOqnYiKIiDuBOyWdGBFeLsiskbKuW7CzCrLugW2vWtPQFyPi28BUSef2Ph4Rn881MjNLZF23YGd53YPCqtY0tCx97apHIGZWQdZ1C3aW1z0orGpNQz9PX+fVLxwzM6u3ak1DPwei0vGI+FguEZmZWV1Vaxq6vG5RmJlZw1RrGnqgnoGYmVljVGsaui0i/lLS42zfRCQgImJC7tGZmVnuqjUNfSF9/Ug9AjGzQaDaeAWPMRiyqjUNrU5fn5e0H3AcyZ3BoxHx+zrFZ2b1Um28gscYDGlZpqH+G+CrJGsPCPi+pH+OiDl5B2dmdVRtvILHGAxpWaahvgR4f0SsA5C0N7CQZKUyMzNrclmmoe4mmWiuxyvAi1kuLmmypKclLZc0o0KZUyQtlrRUknsqmZnVWbVeQxenb18iWY3sTpJnBFOAR/q6sKQW4Grgv5Akk0clzY+IJ0vK7AlcA0yOiBckjRpwTczMbECqNQ2NTF9/l/70uDPjtY8DlvesbibpFpIk8mRJmanAHRHxAkBErMl4bTMzq5FqvYZm7uS1R7N9E1I3cHyvMocBIyTdT5J4vhsRN/S+kKROoBPgoIMO2smwzMysVJZeQ/sAXyRZrWzbymQR8aG+Ti2zr/fcRcOBY4DTgF2B/5D0cEQ8s91JEbOB2QAdHR0V5z8yM7P+y/Kw+CfAU8AYYCawEng0w3ndwIEl223AqjJl7omI1yLiZeBBYGKGa5uZWY1kSQR7R8T1wKaIeCAiPgOckOG8R4GxksZI2gU4B5jfq8ydwMmShkvajaTpaBlmZlY3WcYRbEpfV0s6g+RbfVtfJ0XEZknTgXuBFmBORCyVdEF6/NqIWCbpHmAJsBW4LiKeGEhFzMxsYLIkgm9IehfwD8D3gT2Av89y8Yi4G7i7175re21/B/hOpmjNzKzm+kwEEXFX+nY9cGq+4ZiZWb31+YxA0iGSfi7pZUlrJN0p6ZB6BGdmZvnL8rD4JuA2YD/gAOCnwM15BmVmZvWT5RmBIuLGku0fpw+BzaxIqq1V0Gy8tsJ2qs01tFf69tfphHG3kAwI+yTwv+sQm5kNFtXWKmg2XlthB4ooP1BX0gqSD/6yI4QjoiHPCTo6OqKrq6sRv9rMhoKeu5r7729kFHUnaVFEdJQ7Vm2uoTH5hWRmZoNFlrmGRgD/Hfhguut+4IcRsaniSWZm1jSyPCz+ATCCZN0AgGnpvr/JKygzM6ufLIng2IgonQjuV5IeyysgMzOrryzjCLZIOrRnIx1MtiW/kMzMrJ6y3BH8I0kX0udIehC9F/h0rlGZmeVpMIyJGERjGaomgnTd4YnAWOBwkkTwVES8VYfYzMxqbzCMiRhkYxmqJoKI2CLpYxFxJclU0WZmza2zs/EfwI2+G+klS9PQQkmzgFuB13p2RsRvc4vKzMzqJksiOCl9/eeSfQH0tWaxmZk1gSyJ4Ox0PWEzMxuCKnYflfRRSWuBJZK6JZ1UqayZmTWvauMIvgmcHBEHAGcB36pPSGZmVk/VmoY2R8RTABHxG0kj6xSTmdnQN5CxDO3tcNVVNQ+lWiIYJeniStsRcUXNozEzK4LBMJahRLVE8D+BkVW2zcxsIAbDWIYS1dYjmFnPQMzMrDGyTDpnZmZDmBOBmVnBORGYmRVcxWcEvXoM7cC9hszMhoZqvYbcQ8jMrADca8jMrOD6nHROUivwWWA80NqzPyI+k2NcZmZWJ1keFt8I7AdMAh4A2oBX8gzKzMzqJ0si+JOI+CfgtYiYB5wBHJVvWGZmVi9ZEsGm9PWPkt4HvAs4OLeIzMysrrIsTDNb0ruBfwLmA+9M35uZ2RCQJRH8KCK2kDwfOCTneMzMrM6yNA2tkDRb0mmS1J+LS5os6WlJyyXNqFLuWElbJH2iP9c3M7OdlyURHA78H+BCYKWkWZI+0NdJklqAq4HTgXHAuZLGVSj3L8C9/QnczMxqo89EEBFvRMRtEfFxoB3Yg6SZqC/HAcsj4rmI2AjcAkwpU+7vgNuBNdnDNjOzWsk06ZykP5N0DfBbkkFlf5nhtNHAiyXb3em+0uuOBs4Eru3j93dK6pLUtXbt2iwhm5lZRllGFq8AFgO3AZdExGsZr13ueUL02r4K+FJEbKn2+CEiZgOzATo6Onpfw8zMdkKWXkMTI2LDAK7dDRxYst0GrOpVpgO4JU0C7wE+LGlzRPzbAH6fmZkNQLVpqL8YEd8Gvilph2/hEfH5Pq79KDBW0hjgJeAcYLsVmyNiTMnvmwvc5SRgZlZf1e4IlqWvXQO5cERsljSdpDdQCzAnIpZKuiA9XvW5gJmZ1Yciqje5S3p/RPxnneLpU0dHR3R1DSg3mZkVlqRFEdFR7liWXkNXSHpK0tclja9xbGZm1mBZxhGcCpwCrCWZd+hxSf8j78DMzKw+Mo0jiIjfR8T3gAtIupJ+NdeozMysbvpMBJKOlHSppCeAWcBCkq6gZmY2BGSafRS4GfiLiOg9DsDMzJpc1USQTgj3u4j4bp3iMTOzOqvaNJSuQ7C3pF3qFI+ZmdVZlqah54F/lzQf2DbPUERckVtUZmZWN1kSwar0ZxgwMt9wzMys3vpMBBExsx6BmJlZY2SZhvrX7Dh9NBHxoVwiMjOzusrSNPSPJe9bgbOAzfmEY2Zm9ZalaWhRr13/LinLUpVmZtYEsjQN7VWyOQw4Btgvt4jMzKyusjQNLSJ5RiCSJqEVwGfzDMrMzOonS9PQmL7KmJlZ86o4sljSsZL2K9k+T9Kdkr7Xq7nIzMyaWLUpJn4IbASQ9EHgMuAGYD0wO//QzMysHqo1DbVExB/S958EZkfE7cDtkhbnH5qZmdVDtTuCFkk9ieI04Fclx7I8ZDYzsyZQ7QP9ZuABSS8DbwAPAUj6E5LmITMzGwIqJoKI+Kak+4D9gV9GRM80E8OAv6tHcGZmlr+qTTwR8XCZfc/kF46ZmdVbpsXrzcxs6HIiMDMrOCcCM7OCcyIwMys4JwIzs4JzIjAzKzgnAjOzgnMiMDMrOCcCM7OCcyIwMys4JwIzs4LLNRFImizpaUnLJc0oc/yvJC1JfxZKmphnPGZmtqPcEoGkFuBq4HRgHHCupHG9iq0A/iwiJgBfxyufmZnVXZ53BMcByyPiuYjYCNwCTCktEBELI+L/pZsPA205xmNmZmXkmQhGAy+WbHen+yr5LPCLcgckdUrqktS1du3aGoZoZmZ5JgKV2Rdl9iHpVJJE8KVyxyNidkR0RETHPvvsU8MQzcwsz7WHu4EDS7bbgFW9C0maAFwHnB4R63KMx8zMysjzjuBRYKykMZJ2Ac4B5pcWkHQQcAcwzSufmZk1Rm53BBGxWdJ04F6gBZgTEUslXZAevxb4KrA3cI0kgM0R0ZFXTGZmtiO9vSZ9c+jo6Iiurq5Gh2Fm1lQkLar0Rdsji83MCs6JwMys4JwIzMwKzonAzKzgnAjMzArOicDMrOCcCMzMCs6JwMys4JwIzMwKzonAzKzgnAjMzArOicDMrOCcCMzMCs6JwMys4JwIzMwKzonAzKzgnAjMzArOicDMrOCcCMzMCs6JwMys4JwIzMwKzonAzKzgnAjMzArOicDMrOCcCMzMCs6JwMys4JwIzMwKzonAzKzgnAjMzArOicDMrOCcCMzMCs6JwMys4JwIzMwKzonAzKzgnAjMzAou10QgabKkpyUtlzSjzHFJ+l56fImko/OMx8zMdpRbIpDUAlwNnA6MA86VNK5XsdOBselPJ/CDvOIxM7Py8rwjOA5YHhHPRcRG4BZgSq8yU4AbIvEwsKek/XOMyczMehme47VHAy+WbHcDx2coMxpYXVpIUifJHQPAq5KeHmBM7wFeHuC5g43rMjgNlboMlXqA69LjvZUO5JkIVGZfDKAMETEbmL3TAUldEdGxs9cZDFyXwWmo1GWo1ANclyzybBrqBg4s2W4DVg2gjJmZ5SjPRPAoMFbSGEm7AOcA83uVmQ+cl/YeOgFYHxGre1/IzMzyk1vTUERsljQduBdoAeZExFJJF6THrwXuBj4MLAdeBz6dVzypnW5eGkRcl8FpqNRlqNQDXJc+KWKHJnkzMysQjyw2Mys4JwIzs4IrTCLoa7qLwUzSgZJ+LWmZpKWSvpDu30vSAknPpq/vbnSsWUhqkfSfku5Kt5u1HntK+pmkp9L/Nyc2cV3+Pv3bekLSzZJam6UukuZIWiPpiZJ9FWOX9OX0c+BpSZMaE/WOKtTjO+nf1xJJ/0vSniXHalaPQiSCjNNdDGabgX+IiCOBE4AL0/hnAPdFxFjgvnS7GXwBWFay3az1+C5wT0QcAUwkqVPT1UXSaODzQEdEvI+kc8c5NE9d5gKTe+0rG3v67+YcYHx6zjXp58NgMJcd67EAeF9ETACeAb4Mta9HIRIB2aa7GLQiYnVE/DZ9/wrJB85okjrMS4vNA/5rYyLMTlIbcAZwXcnuZqzHHsAHgesBImJjRPyRJqxLajiwq6ThwG4k43maoi4R8SDwh167K8U+BbglIt6KiBUkPRaPq0ugfShXj4j4ZURsTjcfJhlrBTWuR1ESQaWpLJqOpIOB9wO/AfbtGXeRvo5qXGSZXQV8Edhasq8Z63EIsBb4UdrMdZ2k3WnCukTES8DlwAsk07usj4hf0oR1KVEp9mb+LPgM8Iv0fU3rUZREkGkqi8FO0juB24GLImJDo+PpL0kfAdZExKJGx1IDw4GjgR9ExPuB1xi8TSdVpe3nU4AxwAHA7pL+urFR5aYpPwskfYWkifgnPbvKFBtwPYqSCJp+KgtJI0iSwE8i4o509//tma01fV3TqPgy+lPgY5JWkjTPfUjSj2m+ekDyN9UdEb9Jt39GkhiasS5/DqyIiLURsQm4AziJ5qxLj0qxN91ngaTzgY8AfxVvD/yqaT2KkgiyTHcxaEkSSVv0soi4ouTQfOD89P35wJ31jq0/IuLLEdEWEQeT/D/4VUT8NU1WD4CI+D3woqTD012nAU/ShHUhaRI6QdJu6d/aaSTPoZqxLj0qxT4fOEfSOySNIVkL5ZEGxJeJpMnAl4CPRcTrJYdqW4+IKMQPyVQWzwC/A77S6Hj6GfsHSG77lgCL058PA3uT9Ih4Nn3dq9Gx9qNOpwB3pe+bsh5AO9CV/n/5N+DdTVyXmcBTwBPAjcA7mqUuwM0kzzY2kXxT/my12IGvpJ8DTwOnNzr+PuqxnORZQM+/+2vzqIenmDAzK7iiNA2ZmVkFTgRmZgXnRGBmVnBOBGZmBedEYGZWcE4EZhVI2lvS4vTn95JeSt+/KumaRsdnVivuPmqWgaRLgVcj4vJGx2JWa74jMOsnSaeUrKVwqaR5kn4paaWkj0v6tqTHJd2TTg2CpGMkPSBpkaR7e6Y/MBsMnAjMdt6hJFNrTwF+DPw6Io4C3gDOSJPB94FPRMQxwBzgm40K1qy34Y0OwGwI+EVEbJL0OMmiLvek+x8HDgYOB94HLEim8qGFZCoBs0HBicBs570FEBFbJW2Ktx+8bSX5NyZgaUSc2KgAzapx05BZ/p4G9pF0IiRTiksa3+CYzLZxIjDLWSTLo34C+BdJj5HMInlSY6Mye5u7j5qZFZzvCMzMCs6JwMys4JwIzMwKzonAzKzgnAjMzArOicDMrOCcCMzMCu7/A/yPLOR1fzcvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Get the data\n",
    "inFile1 = 'altman_13_2.txt'\n",
    "inFile2 = 'altman_13_3.txt'\n",
    "url_base = 'https://raw.githubusercontent.com/thomas-haslwanter/statsintro_python/master/ipynb/Data/data_altman/'\n",
    "url1 = url_base + inFile1\n",
    "url2 = url_base + inFile2\n",
    "data_1 = np.genfromtxt(urlopen(url1), delimiter=',')\n",
    "data_2 = np.genfromtxt(urlopen(url2), delimiter=',')\n",
    "\n",
    "# Determine the Kaplan-Meier curves\n",
    "(p1, r1, t1, sp1,se1) = kaplanmeier(data_1)\n",
    "(p2, r2, t2, sp2,se2) = kaplanmeier(data_2)\n",
    "\n",
    "# Make a combined plot for both datasets\n",
    "plt.step(t1,sp1, where='post')\n",
    "plt.step(t2,sp2,'r', where='post')\n",
    "\n",
    "plt.legend(['Data1', 'Data2'])\n",
    "plt.ylim(0,1)\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('Survival Probability')\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X^2 = 3.207\n",
      "p=0.0733, the two survival curves are not signifcantly different.\n"
     ]
    }
   ],
   "source": [
    "'''Logrank hypothesis test, comparing the survival times for two different datasets'''\n",
    "\n",
    "times_1 = data_1[:,0]\n",
    "censored_1 = data_1[:,1]\n",
    "atRisk_1 = np.arange(len(times_1),0,-1)\n",
    "failures_1 = times_1[censored_1==0]\n",
    "\n",
    "times_2 = data_2[:,0]\n",
    "censored_2 = data_2[:,1]\n",
    "atRisk_2 = np.arange(len(times_2),0,-1)\n",
    "failures_2 = times_2[censored_2==0]\n",
    "\n",
    "failures = np.unique(np.hstack((times_1[censored_1==0], times_2[censored_2==0])))\n",
    "num_failures = len(failures)\n",
    "r1 = np.zeros(num_failures)\n",
    "r2 = np.zeros(num_failures)\n",
    "r  = np.zeros(num_failures)\n",
    "f1 = np.zeros(num_failures)\n",
    "f2 = np.zeros(num_failures)\n",
    "f  = np.zeros(num_failures)\n",
    "e1 = np.zeros(num_failures)\n",
    "f1me1 = np.zeros(num_failures)\n",
    "v = np.zeros(num_failures)\n",
    "\n",
    "for ii in range(num_failures):\n",
    "    r1[ii] = np.sum(times_1 >= failures[ii])\n",
    "    r2[ii] = np.sum(times_2 >= failures[ii])\n",
    "    r[ii] = r1[ii] + r2[ii]\n",
    "    \n",
    "    f1[ii] = np.sum(failures_1==failures[ii])\n",
    "    f2[ii] = np.sum(failures_2==failures[ii])\n",
    "    f[ii] = f1[ii] + f2[ii]\n",
    "    \n",
    "    e1[ii] = r1[ii]*f[ii]/r[ii]\n",
    "    f1me1[ii] = f1[ii] - e1[ii]\n",
    "    v[ii] = r1[ii]*r2[ii]*f[ii]*(r[ii]-f[ii]) / ( r[ii]**2 *(r[ii]-1) )\n",
    "\n",
    "    O1 = np.sum(f1)\n",
    "    O2 = np.sum(f2)\n",
    "    E1 = np.sum(e1)\n",
    "    O1mE1 = np.sum(f1me1)\n",
    "    V = sum(v)\n",
    "    \n",
    "chi2 = (O1-E1)**2/V\n",
    "p = stats.chi2.sf(chi2, 1)\n",
    "\n",
    "print('X^2 = {0:5.3f}'.format(chi2))\n",
    "if p < 0.05:\n",
    "    print('p={0:6.4f}, the two survival curves are signifcantly different.'.format(p))\n",
    "else:\n",
    "    print('p={0:6.4f}, the two survival curves are not signifcantly different.'.format(p))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
